#### Physicians in Congress Replication Script 

rm(list = ls()) 
######################################################################
#                   Loading Packages 
######################################################################
library(tidyverse) 
library(ggplot2)
library(readr)
library(ggpubr)
library(sandwich)
library(lmtest)
library(stargazer)
library(plm)
library(haven)
library(DirectEffects)
library(dotwhisker)
library(dplyr)
library(readxl)

setwd("Desktop/Physicians in Congress")
######################################################################
#                   Loading Data
######################################################################
pic <- read_csv("pic.csv") # Load Physician In Congress Data 
pic$doctor <- 1 # physician indicator

cel1 <- read_csv("CELHouse93to116_health_ILES.csv") # LES 93-116
full <- merge(cel1, pic, by="icpsr", all=T) # merge data 
full$doctor[is.na(full$doctor)] <- 0  # code all non-physicians as 0


######################################################################
#               Descriptive  Plots 
######################################################################

full$doctor2 <- NA 
full$doctor2[full$doctor==1] <- "Physician"
full$doctor2[full$doctor==0] <- "Non-Physician"

full$hiles <- log(full$iles_health +1) # logged +1 health LES


ggplot(full, aes(iles_health)) +
  geom_density() + theme_classic() + labs(y="Density", x="Health ILES") # Demonstrates skew of HLES (appendix figure)


# Calculates Average Health LES, Logged Health LES, 
#Health Bill Sponsorship, Health Bills Pass the Chamber for PMCs vs. Others
avgs  <-full %>%
  group_by(doctor2) %>%
  summarise(avgers= mean(iles_health, na.rm=T))

avgs1  <-full %>%
  group_by(doctor2) %>%
  summarise(avgers= mean(hiles, na.rm=T))


avgs2  <-full %>%
  group_by(doctor2) %>%
  summarise(avgers= mean(health_all_bill, na.rm=T))


avgs3  <-full %>%
  group_by(doctor2) %>%
  summarise(avgers= mean(health_all_pass, na.rm=T))


# Plots the averages in outcomes by PMC
p1=ggplot(avgs1, aes(x = doctor2, y =avgers))+
  geom_bar(stat='identity')+
  labs(x="Member", y = "Average Logged H-LES") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_bw()
p1


p2=ggplot(avgs2, aes(x = doctor2, y =avgers))+
  geom_bar(stat='identity')+
  labs(x="Member", y = "Average Health Bills Sponsored") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_bw()
p2

p3=ggplot(avgs3, aes(x = doctor2, y =avgers))+
  geom_bar(stat='identity')+
  labs(x="Member", y = "Average Health Laws Passed") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_bw()
p3


# Combines p1, p2, and p3 to create Figure 2 in Main Manuscript 
ggarrange(p2, p3, p1, ncol=3, nrow=1)


# Calculates # PMCS by Party over time 
trends  <-full %>%
  group_by(congress, Party.x) %>%
  summarise(doc= sum(doctor, na.rm=T))


trends <- trends[which(trends$congress<=116),]
trends <- trends[which(trends$Party.x!="Independent"),]

# Creates Figure 1 in Main Text 
ggplot(trends,aes (x=congress, y=doc, color=Party.x, fill=Party.x)) + xlim(90, 120) + 
  geom_line(size=1)+ 
    scale_fill_grey() + 
    scale_color_grey() + theme_bw() + labs(y="Number of Doctors in Congress", x="Congress", color="Party", fill="Party")




######################################################################
#    Bivariate Regression Results (Health Legislative Effectiveness)
######################################################################

m1 <- lm(health_all_bill ~ doctor, data=full)
m2 <- lm(health_all_pass ~ doctor, data=full)
m3 <- lm(hiles ~ doctor, data=full)


stargazer(m1, m2, m3, 
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


# robust standard errors 
m1robust <- coeftest(m1, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2, 
                     vcov = vcovHC,
                     type = "HC1")



m3robust <- coeftest(m3, 
                     vcov = vcovHC,
                     type = "HC1")


stargazer(m1robust, m2robust, m3robust, 
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 



# Double-checking pre-covid robustness
precovid <- full[which(full$congress<=115),] # Pre-covid (93-115) data 


m1 <- lm(health_all_bill ~ doctor, data=precovid)
m2 <- lm(health_all_pass ~ doctor, data=precovid)
m3 <- lm(hiles ~ doctor, data=precovid)


stargazer(m1, m2, m3, 
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


# robust standard errors 

m1robust <- coeftest(m1, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2, 
                     vcov = vcovHC,
                     type = "HC1")



m3robust <- coeftest(m3, 
                     vcov = vcovHC,
                     type = "HC1")


stargazer(m1robust, m2robust, m3robust, 
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 




######################################################################
#       Full  Regression Results (Health Legislative Effectivess)
######################################################################
# load control variables 
LES <- read_dta("CELHouse93to117ReducedClassic.dta")

full2 <- merge(full, LES, by=c("icpsr", "congress"), all=T) # merge with full data
full2precovid <- merge(precovid, LES, by=c("icpsr", "congress"), all=T) # merge with pre-covid data


# Models Estimated in Table 2 
m1c <- lm(health_all_bill ~ doctor  + dwnom1 + female + majority + power + chair + subchr + seniority +as.factor(congress), data=full2)
m2c <- lm(health_all_pass ~ doctor+ dwnom1 + female + majority + power + chair + subchr + seniority+as.factor(congress), data=full2)
m3c <- lm(hiles~ doctor+ dwnom1 + female + majority + power + chair + subchr + seniority+as.factor(congress), data=full2)


stargazer(m1c, m2c, m3c, 
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Cmte", "Cmte Chair", "Subcmte Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 

# robust standard errors 
m1robust <- coeftest(m1c, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2c, 
                     vcov = vcovHC,
                     type = "HC1")



m3robust <- coeftest(m3c, 
                     vcov = vcovHC,
                     type = "HC1")


# results presented in Table 2 
stargazer(m1robust, m2robust, m3robust,
          dep.var.labels=c("Sponsor Health Bill", "Pass Health Bill", "Logged Health LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Cmte", "Cmte Chair", "Subcmte Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 



# Robusness check with raw health LES 
m4c <- lm(iles_health ~ doctor + dwnom1 + female + majority + power + chair + subchr + seniority+as.factor(congress), data=full2) # (not in article, for reviewer)


m4robust <- coeftest(m4c, 
                     vcov = vcovHC,
                     type = "HC1")

stargazer(m4c,
          dep.var.labels=c("Raw Health LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Cmte", "Cmte Chair", "Subcmte Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 

stargazer(m4robust,
          dep.var.labels=c("Raw Health LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Cmte", "Cmte Chair", "Subcmte Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


# Probing substantive significance of results (majority party vs PMC status)

majority <- full2[which(full2$majority==1),]


# PMC Average versus majority party average 
difs1  <-majority %>%
  dplyr::group_by(congress, doctor) %>%
  dplyr::summarise(meanloghiles= mean(hiles, na.rm=T))

difs2  <-majority %>%
  dplyr::group_by(doctor) %>%
  dplyr::summarise(meanloghiles= mean(hiles, na.rm=T))


difs1 <- difs1[which(difs1$congress<=116),]

difs1 <- difs1[complete.cases(difs1), ]


difs1$PMC <- NA
difs1$PMC[difs1$doctor==1] <- "PMC"
difs1$PMC[difs1$doctor==0] <- "Majority"

# Creates Figure 10 in Appendix showing PMC verses majority party average logged health LES 
p1 <- ggplot(difs1, aes(congress, meanloghiles, color=PMC, fill=PMC)) +
  geom_bar(stat="identity", position="dodge") + scale_color_grey() +
  scale_fill_grey() + labs(y="Average Logged Health LES Score", x= "Congress") + 
  theme_classic() + annotate(geom="text", x=110, y=2.5, label="PMC Average = 0.37",
                             color="grey", size=3)+ annotate(geom="text", x=110, y=2.25, label="Majority Average = 0.32",
                                                             color="black", size=3)

p1




######################################################################
#       Full Regression Results (Overall Legislative Effectiveness)
######################################################################

# Models used in Table 6 (appendix) used to create plotted "overall" results for Figure 6 in Main Text

m1d <- lm(allbill1 ~ doctor  + dwnom1 + female + majority + power + chair + subchr + seniority +as.factor(congress), data=full2)
m2d <- lm(allpass1 ~ doctor+ dwnom1 + female + majority + power + chair + subchr + seniority+as.factor(congress), data=full2)
full2$logles <- log(full2$les + 1)
m3d <- lm(logles ~ doctor+ dwnom1 + female + majority + power + chair + subchr + seniority+as.factor(congress), data=full2)

stargazer(m1d, m2d, m3d, 
          dep.var.labels=c("Sponsor Bill", "Pass Bill", "LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Cmte", "Cmte Chair", "Subcmte Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 

# robust standard errors 
m1robust <- coeftest(m1d, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2d, 
                     vcov = vcovHC,
                     type = "HC1")



m3robust <- coeftest(m3d, 
                     vcov = vcovHC,
                     type = "HC1")


stargazer(m1robust, m2robust, m3robust,
          dep.var.labels=c("Sponsor Bill", "Pass Bill", "Logged LES"),
          covariate.labels=c("Physician", "DW Nominate (Dim1)", "Female",
                             "Majority",  "Power Committee", "Committee Chair", "Subcommittee Chair", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


######################################################################
#       Full Regression Results (Ideology Analyses)
######################################################################


full2$dems <- NA 
full2$dems[full2$Party.x=="Democrat"]<- 1 
full2$dems[full2$Party.x=="Republican"]<- 0 # Party indicator 

dems <- full2[which(full2$Party.x=="Democrat"),]
gops <- full2[which(full2$Party.x=="Republican"),] # split data by party 

# Models used to create Figure 3 in main text and Table 7 in the appendix 
m1 <- lm(dwnom1 ~ doctor + dems + female + seniority +as.factor(congress), data=full2)
m2 <- lm(dwnom1 ~ doctor + female + seniority +as.factor(congress), data=dems)
m3 <- lm(dwnom1 ~ doctor + female + seniority +as.factor(congress), data=gops)

# robust standard errors 
m1robust <- coeftest(m1, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2, 
                     vcov = vcovHC,
                     type = "HC1")

m3robust <- coeftest(m3, 
                     vcov = vcovHC,
                     type = "HC1")


# Table 7 
stargazer(m1robust,  m2robust, m3robust,
          dep.var.labels=c("DW Nominate"),
          covariate.labels=c("Physician", "Party", "Female", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


stargazer(m1, m2, m3, 
          dep.var.labels=c("DW Nominate"),
          covariate.labels=c("Physician", "Party", "Female", "Seniority"), style = "ajps", 
          omit = c("congress"),
          omit.stat=c("LL","ser","f"), no.space=FALSE, align=TRUE) 


# Figure 3 
dv1 <- dwplot(list(m1robust,  m2robust, m3robust),
              vline = geom_vline(
                xintercept = 0,
                colour = "grey60",
                linetype = 2
              ),
              vars_order = c("doctor"),
              model_order = c("Model 3", "Model 2", "Model 1")
) %>% # plot line at zero _behind_coefs
  relabel_predictors(
    c(doctor = "Physician")
  ) +
  theme_bw(base_size = 4) + 
  # Setting `base_size` for fit the theme
  # No need to set `base_size` in most usage
  xlab("DW-NOMINATE") + ylab("") +
  geom_vline(xintercept = 0,
             colour = "grey60",
             linetype = 2) +
  theme(
    legend.position = "bottom", 
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_blank()
  ) +
  scale_colour_grey(
    start = .8,
    end = .4,
    name = "Model",
    breaks = c("Model 1", "Model 2", "Model 3"),
    labels = c("All", "Democrats", "Republicans")
  ) +
  scale_fill_grey(
    start = .8,
    end = .4,
    name = "Model",
    breaks = c("Model 1", "Model 2", "Model 3"),
    labels = c("All", "Democrats", "Republicans")
  ) +theme(axis.text=element_text(size=12),
           axis.title=element_text(size=14,face="bold"))+
  theme(legend.key.height= unit(1, 'cm'),
        legend.key.width= unit(1, 'cm'),
        legend.title=element_text(size=10), 
        legend.text=element_text(size=9))


######################################################################
#       Full Regression Results (Health Committees Appendix)
######################################################################

cmtemoney <- read_dta("How_Do_interest_Groups_Seek_Access_to_Committees.dta") # health committee data 

full2$name <- toupper(full2$thomas_name.x)

full3 <- merge(full2, cmtemoney, by.x = "name",  by.y = "Klarner_name") # merge with full data 

# limit data to committee range 

full3 <- full3[which(full3$congress<=113 & full3$congress>=100),]

summary(full)
# Appendix materials 
# Probability of Health committee assignment 
c1 <- lm(cmt_health ~ doctor + dwnom1 + female + power + chair + subchr +  seniority +as.factor(congress), data=full3)

# Donations from health organizations 
c2<- lm(log_amount36_1 ~ doctor*cmt_health + dwnom1 + female + majority  + power + chair + subchr +  seniority + as.factor(congress), data=full3)
c3<- lm(log_amount37_1 ~ doctor*cmt_health + dwnom1 + female + majority  + power + chair + subchr +  seniority+ as.factor(congress), data=full3)
c4<- lm(log_amount39_1 ~ doctor*cmt_health+ dwnom1 + female + majority  + power + chair + subchr +  seniority+ as.factor(congress), data=full3)


# robust standard errors 
c1robust <- coeftest(c1, 
                     vcov = vcovHC,
                     type = "HC1")

c2robust <- coeftest(c2, 
                     vcov = vcovHC,
                     type = "HC1")

c3robust <- coeftest(c3, 
                     vcov = vcovHC,
                     type = "HC1")

c4robust <- coeftest(c4, 
                     vcov = vcovHC,
                     type = "HC1")



# Health Policy Effectiveness of physicians on health committees. 

healths<- full3[which(full3$cmt_health==1),]

m1da <- lm(health_all_bill ~ doctor  + dwnom1 + female + majority  + power + chair + subchr + seniority, data=healths)
m2da <- lm(health_all_pass ~ doctor+ dwnom1 + female + majority + power + chair + subchr + seniority, data=healths)
m4da <- lm(hiles~ doctor +  dwnom1 + female + majority  +  power + chair + subchr + seniority, data=healths)


m1robust <- coeftest(m1da, 
                     vcov = vcovHC,
                     type = "HC1")


m2robust <- coeftest(m2da, 
                     vcov = vcovHC,
                     type = "HC1")

m3robust <- coeftest(m4da, 
                     vcov = vcovHC,
                     type = "HC1")



######################################################################
#                 Bill Content Analyses Plots 
######################################################################

sponsors <- read_csv("bill_sponsor_categories.csv")

sponsors3 <- sponsors[which(sponsors$Party!="Overall"),]

# Saved Calculations from hand coding for Figure 5
p=ggplot(sponsors3, aes(x = Policy, y =Percent, fill=Party))+
  geom_bar(stat='identity', position="dodge")+
  labs(x="", y = "Percentage", fill="Party") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
p

# Saved Calculations from hand coding  for  for Figure 4
codes_comparison <- read_csv("codes_comparison.csv")
p=ggplot(codes_comparison, aes(x = Category, y =Percent, fill=Group))+
  geom_bar(stat='identity', position="dodge")+ ylim(0, 35)+ 
  labs(x="", y = "Percentage", fill="Group") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
p





######################################################################
#                   Remaining Plots and Tables 
######################################################################

# Creates Descriptive Statistics Table 1 (Main Text abridged) and  Descriptive Statistics Table 3 (appendix full)
full3 <- full2[which(full2$congress<=116),]

full3$logles <- log(full3$les +1)

keep <- c("doctor", "allbill1", "allpass1", "logles", "health_all_bill", "health_all_pass", "hiles", "dwnom1", "female", "majority", "power", "chair", "subchr", "seniority")

full3 <- full3[keep]
stargazer(full3)



# Health committee probability by PMC status 

full3$doctor2 <- NA 
full3$doctor2[full3$doctor==1] <- "Physician"
full3$doctor2[full3$doctor==0] <- "Non-Physician"



results  <-full3 %>%
  group_by(doctor2) %>%
  summarise(avgers= mean(cmt_health, na.rm=T))


# Creates Figure 7 
p1=ggplot(results, aes(x = doctor2, y =avgers))+
  geom_bar(stat='identity')+
  labs(x="Physician Status", y = "Pr(Health Committee Assignment)") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + theme_classic() + ylim(0,0.4)
p1




# Creates Figure 6 Based On Saved Regression Results From Above 

effects <- read_csv("HEALTH_REGLES.csv")

effects$Effect
limits = aes(ymax = Effect + (1.96*SD), ymin=Effect - (1.96*SD))

dodge = position_dodge(width=0.9)

apatheme=theme_bw()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Times'))



p=ggplot(effects, aes(x = Outcome, y = Effect, color=Group, fill=Group))+
  geom_bar(stat='identity', position="dodge")+
  geom_errorbar(limits, position=dodge, width=0.25, color="black")+
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Effect of Physician Status', color="Type of Policy", fill="Type of Policy") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p


effects2 <- read_csv("bills_predicts.csv")


limits = aes(ymax = Effect + (1.96*SD), ymin=Effect - (1.96*SD))

dodge = position_dodge(width=0.9)

apatheme=theme_bw()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        axis.line=element_line(),
        text=element_text(family='Times'))


p1=ggplot(effects2, aes(x = Outcome, y = Effect, color=Group, fill=Group))+
  geom_bar(stat='identity', position="dodge")+
  geom_errorbar(limits, position=dodge, width=0.25, color="black")+
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Predicted Sponsorships', color="Member", fill="Member", x="Type of Bill Sponsored") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p1

# Assembles combined Figure 6
ggarrange(p, p1, ncol=1, nrow=2)



# Generates subfigures for Figure 8  based on saved results 
effects <- read_csv("DIFFERENCE_BYCOMMITTEE.csv")

p1=ggplot(effects, aes(x = Committee, y = HealthBills, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ ylim(0, 15) + 
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Health Sponsor', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p1

p2=ggplot(effects, aes(x = Committee, y = HealthPassed, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ylim(0, 2) + 
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Health Passed', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p2

p3=ggplot(effects, aes(x = Committee, y = HealthLES, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ ylim(0, 1) + 
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Logged Health LES', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p3

p4=ggplot(effects, aes(x = Committee, y = AllBills, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ylim(0, 15) +
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Bills Sponsored', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p4

p5=ggplot(effects, aes(x = Committee, y = AllPassed, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ylim(0, 2) + 
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Bills Passed', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p5


p6=ggplot(effects, aes(x = Committee, y = LLES, color=Doctor, fill=Doctor))+
  geom_bar(stat='identity', position="dodge")+ ylim(0, 1) + 
  apatheme + geom_hline(yintercept = 0, linetype="dashed", color="black") + 
  labs(y='Logged LES', color="Member", fill="Member", x="Health Committee Assignment") + theme(text = element_text(size = 17)) + 
  scale_fill_grey() + scale_color_grey() 
p6

# Assembles  Figure 8 based on saved results 
ggarrange(p1,  p4, p2, p5, p3, p6, common.legend = T, ncol = 2, nrow = 3, legend = "bottom")


